This document represents the Pytho codes associated with
Series Representations and Simulations of Isotropic Random Fields in the Euclidean Space
by Zhengwei Ma and Chunsheng Ma
# python imports
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import itertools
from numpy.random import normal, uniform, gamma
from scipy.special import jv
import plotly.graph_objects as go
import multiprocessing as mp
Here we set some initial values for simulation, as well as create some helper functions
r = 10
theta = 0.2
simulation_number = 1000 # number of times to simulate
m = 100 # number of points to simulate
r_list = np.linspace(0, 10, m)
theta_list = np.linspace(0, np.pi * 2, m)
R, T = np.meshgrid(r_list, theta_list) # R, T are matrices of size 100 by 100
X, Y = R*np.cos(T), R*np.sin(T) # express in polar coordinates
def generate_z(function):
z_list = []
with mp.Pool(mp.cpu_count()) as p:
z_list = p.starmap(function, itertools.product(range(m), range(m)))
Z = np.reshape(z_list, (m, m))
return Z
We want \begin{equation} y_{1} = -log(U_{0, 1}) \end{equation}
\begin{equation} y_{2} = U_{log(2), log(8)} \end{equation}\begin{equation} v_{1} = y_{1}^{0.5} * exp(\frac{-y_{2}}{2}) \end{equation}\begin{equation} w_{1} = U_{0, 1} \end{equation}\begin{equation} u_{1} = U_{0, 1} \end{equation}def example_1(i, j):
"""
This function takes in two positional arguments, i and j, which represents the corresponding
values in a m by m matrix (based on the value set in cell 4)
"""
r = R[i, j]
theta = T[i, j]
z_output = []
for i in range(simulation_number):
y1 = -np.log(uniform(low=0, high=1, size=1))
y2 = uniform(low=np.log(2), high=np.log(8), size=1)
v1 = (y1 ** 0.5) * np.exp(-y2/2)
w1 = uniform(low=0, high=1, size=1)
u1 = uniform(low=0, high=1, size=1)
Zi = jv(i, 2*r*np.sqrt(-v1 * np.log(w1))) * np.cos(i*theta+2*np.pi*u1)
z_output.append(Zi)
return np.sqrt(2) * np.sum(z_output)
Z = generate_z(example_1)
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z=dict(
show=True, usecolormap=True,
highlightcolor="limegreen"))
fig.update_layout(title='', autosize=False,
width=1000, height=1000,
margin=dict(l=65, r=50, b=65, t=90), scene=dict(
xaxis=dict(showticklabels=False),
yaxis=dict(showticklabels=False),
zaxis=dict(showticklabels=False),
))
fig.show()
We want \begin{equation} y_{1} = \Gamma(0.5, 1) \end{equation}
\begin{equation} y_{2} = U_{2^{0.5}, 8^{0.5}} \end{equation}\begin{equation} v_{1} = \frac{y_{1}^{0.5}}{y_{2}} \end{equation}\begin{equation} w_{1} = U_{0, 1} \end{equation}\begin{equation} u_{1} = U_{0, 1} \end{equation}def example_2(i, j):
r = R[i, j]
theta = T[i, j]
z_output = []
for i in range(simulation_number):
y1 = gamma(shape=0.5, scale=1, size=1)
y2 = uniform(low=(2**0.5), high=(8**0.5), size=1)
v1 = (y1 ** 0.5) / y2
w1 = uniform(low=0, high=1, size=1)
u1 = uniform(low=0, high=1, size=1)
Zi = jv(i, 2*r*np.sqrt(-v1 * np.log(w1))) * np.cos(i*theta+2*np.pi*u1)
z_output.append(Zi)
return np.sqrt(2) * np.sum(z_output)
Z = generate_z(example_2)
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z=dict(
show=True, usecolormap=True,
highlightcolor="limegreen"))
fig.update_layout(title='', autosize=False,
width=1000, height=1000,
margin=dict(l=65, r=50, b=65, t=90), scene=dict(
xaxis=dict(showticklabels=False),
yaxis=dict(showticklabels=False),
zaxis=dict(showticklabels=False),
))
fig.show()
We want \begin{equation} y_{1} = \Gamma(3, 1) \end{equation}
\begin{equation} y_{2} = U_{\frac{1}{64}, \frac{1}{4}} \end{equation}\begin{equation} v_{1} = y_{1}^{0.5} * y_{2}^{0.25} \end{equation}\begin{equation} w_{1} = U_{0, 1} \end{equation}\begin{equation} u_{1} = U_{0, 1} \end{equation}def example_3(i, j):
r = R[i, j]
theta = T[i, j]
z_output = []
for i in range(simulation_number):
y1 = gamma(shape=3, scale=1, size=1)
y2 = uniform(low=(1/64), high=(1/4), size=1)
v1 = (y1 ** 0.5) * (y2 ** 0.25)
w1 = uniform(low=0, high=1, size=1)
u1 = uniform(low=0, high=1, size=1)
Zi = jv(i, 2*r*np.sqrt(-v1 * np.log(w1))) * np.cos(i*theta+2*np.pi*u1)
z_output.append(Zi)
return np.sqrt(2) * np.sum(z_output)
Z = generate_z(example_3)
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z=dict(
show=True, usecolormap=True,
highlightcolor="limegreen"))
fig.update_layout(title='', autosize=False,
width=1000, height=1000,
margin=dict(l=65, r=50, b=65, t=90), scene=dict(
xaxis=dict(showticklabels=False),
yaxis=dict(showticklabels=False),
zaxis=dict(showticklabels=False),
))
fig.show()
def besselnum(n):
results = []
i = 0
while i < n:
u = uniform(0, 1)
ug = uniform(0, 1)
g = 2 * ug/(1-ug)
if 3 * (0.5 / ((1+0.5*g) ** 2) * u) <= 2/g*(jv(g, 1) ** 2):
i+=1
results.append(g)
return results
def example_4(i, j):
r = R[i, j]
theta = T[i, j]
u = uniform(low=0, high=1, size=simulation_number)
v = besselnum(simulation_number)
z_list = []
for i in range(simulation_number):
Z = jv(i, r * v[i]) * np.cos(i*theta+2*np.pi*u[i])
z_list.append(Z)
return np.sqrt(2)*sum(z_list)
Z = generate_z(example_4)
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z=dict(
show=True, usecolormap=True,
highlightcolor="limegreen"))
fig.update_layout(title='', autosize=False,
width=1000, height=1000,
margin=dict(l=65, r=50, b=65, t=90), scene=dict(
xaxis=dict(showticklabels=False),
yaxis=dict(showticklabels=False),
zaxis=dict(showticklabels=False),
))
fig.show()